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ABSTRACT 


In this thesis a stochastic model for the data acquisition system in a 
multispectral scanner system, like the one utilized by the LANDSAT satellites, 
is presented. A list of noise sources which are known or presumed to have a 
significant effect in the information extraction process was constructed. Since 
the shot noise introduced by the photodetectors in the sensor system is signal 
level dependent, an atmospheric model was adopted which could adequately 
describe the amount of radiation that gets into the sensors based on the 
atmospheric transmittance. An analysis was carried out to find the output 
spectral statistics in terms of the input signal statistics and the system 
parameters. This was integrated into a set of fortran programs that when 
supplied with, the class statistics, the noise levels introduced by the sensor 
system, the atmospheric transmittance, and the atmospheric path radiance, can 
be used to estimate the classification performance. In order to show the 
beneficts of this model a series of runs were performed in which the Thematic 
Mapper multispectral scanner was the system under consideration. 
Consideration was given to the usage of preprocessing spatial filters as a way to 
combat the noise introduced into the signal at different stages of the system. 
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CHAPTER 1 
INTRODUCTION 


1.1 Background 

Remoted Sensing of the environment is largely concerned with the 
measurements of electromagnetic energy emitted or reflected from the earth 
surface. The data acquisition system can be considered in four basic parts: the 
radiation source, the atmospheric path, the target, and the sensor. In the case 
of Landsat, the sun is the radiation source and the atmosphere modifies the 
spectral distribution of the solar radiance. A portion of the spectral radiance 
that arrives at the earth surface in the visible and reflective infrared is reflected 
back through the atmosphere and then may be measured in several spectral 
bands by a multispectral sensor. Some of the sun’s energy is absorbed and 
then re-emitted at thermal infrared wavelengths. The output of the sensor in 
each spectral band for a pixel on the earth surface is used to form a point in a 
n dimensional space. This data is digitized and transmitted to a ground 
station for processing. 

Different ground covers have different spectral reflectance properties. This 
provides the basis for their identification. Detecting these spectral differences 
between ground covers allows classification of each pixel of the observed scene 
as coming from one of a set of possible classes. Due to the randomness of 
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nature the points from a particular class can be characterized by a stochastic 
ensemble. Fu, Landgrebe, and Phillips [F3] have shown empirically that a 
multivariate gaussian density function is a good characterization for 
agricultural multispectral data. This has imposed the use of statistical 
pattern-recognition methods in the data analysis. A commonly used pattern 
classifier algorithm is the minimum Bayes error scheme [S3], 


1.2 Statement of the problem 

The signal in a remote sensing system is corrupted, at different stages of 
the data acquisition process, by different noise sources. It would be important 
to understand the manner in which the different types of noise occurring in a 
remote sensing system impact the performance of the data analysis portion of 
the system. This subject has received relatively little attention to this time, 
and the amount of information available in the literature is rather limited. 
Studies to date have tended to concentrate primarily on the effect of 
independent white Gaussian noise on classification accuracy [Wl,Ml,M2]. 
Mobasseri developed a parametric model to analytically evaluate the response 
of a multispectral scanner in any operational environment. In his work the 
atmospheric effect was not considered, and the noise introduced by the sensor 
system was arbitrarily chosen; it was assumed to have the same statistical 
properties over all the spectral bands. Maxwell had also made a system analysis 
study of the remote sensing system. He discussed the usage of preprocessing 
the data before any classification is done. Although there have been extensive 
studies on models for ’’correcting” or calibrating data in the face of 
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atmospheric effects [Tl], little is known about the deleterious impact of these 
effects on the process of information extraction. Other system noise sources 
have received even less attention; thus the relatively broad field of noise 
sources, the interrelation of them with remote sensing system parameters, and 
the information extraction process remains unexplored. 

If the full potential of current, and especially of planned sensor systems is 
to be realized, a better delineation of the relationship between sensor 
parameters and data analysis results is needed. 

1.3 Objectives 

It is the objective of this work to develop a stochastic system model that 
can be used to integrate and investigate the effect that different noise sources, 
introduced by the atmosphere and the sensor system, can have on system 
performance. Since both the signal and the noise are considered to be random 
in character and therefore not obviously identifiable by inspection, a suitable 
definition of the two is needed. The signal is defined to be that portion of the 
scene response variability which contributes in any way to the ability to 
discriminate among classes. For example, texture of a given cover class in a 
given scene may appear to be random. However, in this case the randomness 
implies variability which might be indicative of the cover class, and thus might 
be part of the signal. All scene and system response variability which does not 
so contribute to identifi ability is defined to be the noise. The theory of 
stochastic process models for such situations provides a rich background of well 
established tools, and these models have been the traditional types of models 
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used for highly complex problems. 

Since the ultimate goal of most systems designed to recognize patterns is 
to do so accurately, the probability of error is an important index of 
performance which it is desirable to minimize, although, there is a lower bound 
on the probability of error imposed by the random nature of the signals. The 
classifier used is a Bayes minimum error pixel classifier for multivariate 
gaussian density functions [S3]. In order to provide a measure of specificity to 
this otherwise very general study, it was decided to use the Thematic Mapper 
sensor parameters as an example system. In Appendix A a general description 
of the Thematic Mapper is presented. The next chapter will discuss the 
different noise sources to be modeled in the system. 
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CHAPTER 2 
NOISE SOURCES 


2.1 Noise List 

As was mentioned earlier we wish to model scene and system response 
variabilities in the data acquisition system which can potentially improve the 
information extraction process. To begin the study a list of noise sources was 
constructed. This list is intended to include all the sources of noise which are 
known or presumed to have a significant deleterious effect on the classification 
of the data. Since it would be impractical to deal with a very extensive list, it 
was decided to prioritize the list, which is as follows: 

Priority 1 

1. Atmospheric effect(constant atmosphere) 

2. Detector noise process 

a. Shot noise(signal level dependent) 

b. Johnson noise 

3. Quantization effect 
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Priority 2 

1. Training statistics estimation error 

2. Optically induced noise 

3. Atmospheric effect (non-constant atmosphere) 

4. Goniometric variations 

5. Non-class conditional earth surface variability 

In constructing this list consideration was given both to the speculated 
magnitude of the noise effect on net system performance, and upon the 
probability of affecting performance improvements by better parameter 
selection in the system design phase. It was shown by Mobasseri, McGillem, 
and Anuta [M2] that the noise prior to detection is quite negligible and thus for 
all practical purposes can be neglected. This implies that random noise 
generated by the detector and quantizer stages are the major sources of 
disturbance on the signal. This work is limited to the Priority 1 list. We 
proceed next with a description of the noise sources in the Priority 1 list. 


2.2 Atmospheric Effect 

In this section the atmospheric model that was used in the simulation is 
presented. The atmosphere is known to be a quite complex portion of the 
system, with a number of mechanisms operative which will affect radiances 
passing through it. To attempt to model all of these would be prohibitively 
complex and indeed is not necessary for the purposes of this early stage of 
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studying the effect of such noise sources. Instead we will use a simple 
atmospheric model for the preliminary work. 

The atmosphere modifies the spectral solar radiance in a wavelength 
dependent manner by two basic mechanisms, absorption and scattering. The 
combination of these two is frequently refered to as attenuation. The 
fundamentals of atmospheric scattering, and absorption involve complexities in 
mathematical physics beyond the scope of this work . The discussion that 
follows will deal with the topic in a less fundamental manner. 

Absorption is the transformation of radiant energy into heat. In a clear 
(haze free) atmosphere, absorption is almost negligible in many portions of the 
optical range. For a hazy or polluted atmosphere absorption plays an 
important role. Atmospheric absorption due to ozone is very strong for 
wavelengths below ,29um . 

Scattering occurs when the radiation is reflected or refracted, from its 
original directional flow, by particles in the atmosphere. There are three basic 
kinds of scattering mechanism: Rayleigh, Mie, and non-selective scattering, 
Each of these is based upon the presence of different scattering elements. 
Certain spectral bands of the solar irradiance that pass through the atmosphere 
are not severely attenuated. These bands are called atmospheric windows and 
are widely used by earth observational satellites. 

The atmosphere has three basic effects on the data acquisition system. 
First, the spectral and spatial distribution of the incoming solar radiance (the 
source) is modified as it passes through the atmosphere. Second, the reflected 
and/or emitted spectral radiance from the earth surface (the target) is 
attenuated. Third, a component of scattered radiation called path radiance is 
added to the attenuated signal. The path radiance also depends on the 
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reflectivity of the target and its surroundings. For simplicity we will assume 
that the fields of interest are big enough such that the reflectivity is fairly 
uniform around the target, even though this assumption is frequently not the 
case. Kiang [Kl] showed that atmospheric effects on the reflected radiance can 
be approximated with a linear equation, where the spectral radiance that falls 
within the instantaneous field of view (IFOV) of the multispectral scanner is as 
a function of wavelength X, is 

L(X) = r a (X)L s (X) + L p (X). (2.1) 


L s (the signal) is the reflected and/or emitted spectral radiance from the pixel 
under observation, r a (X) is the spectral atmospheric transmittance, and L p is 
the path spectral radiance introduced by scattering. Both L s and L are known 
to be random variables that are functions of the wavelength. Although it is 
not written explicitly r a , L s , and L p are functions of the solar zenith angle (0) 
and the metereological range or visibility (V^). We are assuming that the 
ground is a Lambertian surface, which implies that the reflected radiance is 
independent of the viewing angle. 

The spectral transmittance from the earth surface into the outer 
atmosphere is given by 




~4.(*.V,)sec0 
C ) 


(2.2) 


where T ext is the extinction optical thickness introduced by ozone absorption, 
Rayleigh, and aerosol scattering. Elterman [El] has done experimental and 
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Table 2.1 Atmospheric extinction optical thickness for eight 
metereological ranges and twenty different wave- 
lengths in the range of .27-2.17 urn. 


Extinction Optical Thickness T. 


ext 


Metereological Range VJkm) 


A lr“ u / 

2 

3 

4 

5 

6 

8 

10 

13 

0.27 

76.276 

75.324 

74.825 

74.513 

74.299 

74.020 

73.847 

73.681 

0.28 

40.658 

39.760 

39.289 

38.995 

38.794 

38.530 

38.368 

38.211 

0,30 

7.657 

6.809 

6.364 

6.086 

5.896 

5.647 

5.493 

5.345 

0.32 

4.080 

3.281 

2.863 

2.601 

2.422 

2.187 

2.042 

1.903 

0.34 

3.414 

2.665 

2.273 

2.027 

1.859 

1.639 

1.503 

1.372 

0.36 

3.086 

2.383 

2.014 

1.783 

1.625 

1.418 

1.291 

1.167 

0.38 

2.881 

2.202 

1.847 

1.624 

1.472 

1.272 

1.149 

1.030 

0.40 

2.593 

1.969 

1.642 

1.437 

1.297 

1.114 

1.000 

.891 

0.45 

2.203 

1.650 

1.360 

1.178 

1.054 

.891 

.791 

.694 

0.50 

1.968 

1.462 

1.197 

1.031 

.917 

.768 

.676 

.587 

0.55 

1.805 

1.337 

1.092 

.939 

.834 

.696 

.611 

.529 

0.60 

1.624 

1.204 

.984 

.846 

.751 

.627 

.551 

.477 

0.65 

1.452 

1.069 

.868 

.742 

.655 

.542 

.473 

.405 

0.70 

1.341 

.982 

.793 

.675 

.594 

.488 

.423 

.359 

0.80 

1.178 

.859 

.692 

.588 

.516 

.422 

.364 

.307 

0.90 

1.067 

.777 

.624 

.529 

.463 

.378 

.325 

.273 

1.06 

.961 

.699 

.561 

.475 

.416 

.338 

.290 

.244 

1.26 

.876 

.637 

.511 

.433 

.379 

.308 

.264 

.222 

1.67 

.753 

.548 

.441 

.373 

.328 

.268 

.228 

.192 

2.17 

.672 

.489 

.392 

.332 

.290 

.235 

.202 

.169 
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Figure 2.1 Atmospheric optical thickness as a function of 
wavelenght and metereological range. 




Figure 2.2 Atmospheric transmittance as a function of mete- 
reological range and wavelength. 
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theoretical tabulations of spectral optical parameters for a hazy atmosphere 
with different metereological ranges. A selected part of Elterman's tabulations 
for T^ xt are presented in Table 2.1. Using this table a satisfactory linear 
interpolation can be made between the wavelengths of .27um and about lum, 
because r^ xt is a slowly varying function of the wavelength. Fig. 2.1 is the 
result of interpolating over the wavelength range of .3-1.24 microns and over 
the metereological range of 2 to 13 kilometers. Fig. 2.2 is the atmospheric 
transmittance obtained from (2.2) for a solar zenith angle of 37.5 degrees and 
optical thickness given by Fig. 2.1. Observe that atmospheric attenuation is 
stronger in the blue wavelength region compared to the red wavelength region. 
This accounts for having reddish sunsets. 

There are different models to estimate the atmospheric path radiance. I 
will adopt in this work Duntley’s method [D2], because of its simplicity 
compared to other existing methods (S2|. Duntley argued that each segment of 
the atmosphere has an equilibrium radiance L eq with a path radiance given by 

L p (X) = L eq (X)|l - r„(X,V„,0)]. (2.3) 


In this method the equilibrium radiance can be determined from ground based 
radiance and transmittance measurements. 

Landsat photodetectors respond to the radiation that lies within certain 
narrow spectral bands. Since we are only using a finite number of spectral 
bands a matrix formulation of (2.1) is required. 


12 


L, 




K 


L p, 

l 2 

— 

© 

o 


L s, 

+ 

Lpa 



T a. 


K 


Lp. 


(2.4) 


L = r a L s + L p . 


Lj is the radiance that enters the pupil of the sensor in the wavelength range 
covered by the i th channel of the sensor. For notational simplicity r a , L p , and 
L s will be refered to from now on as the transmittance matrix, the equivalent 
radiance vector and the signal vector, respectively. 

A linear transform of the form of (2.4) does not of itself affect the 
classification results obtained from a minimum error Bayes classifier for 
multivariate gaussian distributions. This can be shown as follows. Denote the 
density function of the reflected radiance for class Wj before it passes through 
the atmosphere as 

P(LJ wj) = N(Mj,£j) (2.5) 

that is, a multivariate gaussian density function with covariance matrix Ej and 
mean vector Mj . The density function of the received radiance is given by 

P(L| wj) = N(r a Mj + L p , r a Ejr‘). (2.6) 


Using (2.4) it is a matter of substitution to show 
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p(L| wj) 


1 



P(L S | w j)- 


(2.7) 


Since the transformed density function differs from its original density by a 
multiplicative constant, the classification statistics in a multiclass problem are 
unaffected. It can also be shown that the correlation coefficients between 
channels is also preserved by (2.4), 

Pi) = Piy (2-8) 

Thus the probability of misclassification between two different classes remain 
the same after going through the atmosphere. This is true as long as an ideal 
sensor and a constant atmosphere are considered. 

A plausible atmospheric model has been introduced in this section. The 
inputs to the model are: the atmospheric transmittance matrix and the path 
radiance vector. It is important to recognize that this model is deterministic in 
nature and does not account for spatial, spectral, and temporal atmospheric 
variations. In Appendix B a non constant atmospheric model is discussed. 
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2.3 Shot Noise 

The purpose of this section is to find the first order statistics of the shot 
noise introduced by the photodetectors in the sensor system. Shot noise 
processes were first studied in connection with phototubes, but they have also 
been observed in solid states devices like photodiodes [Dlj. 

A photodetector converts light intensity into a current or voltage 
waveform. When light with the proper wavelength is received at the surface of 
the photosensitive material, electrons are released from its inner surface and 
pulled by an electric field to a collecting anode. The flow of a single electron 
from the cathode to the anode produces a current pulse, h(t), at the 
photodetector output. This pulse is of finite duration and area equal to one 
electron charge [Dl]. If an electromagnetic field is received at the surface of 
the photodetector at time t^O, the output response current of the 
photodetector is the sum of a random number of randomly located response 
functions h(t), 

k(0,t) 

x(t) = g h(t-t n ) t>0, (2.9) 

n=0 


where t n is the time of release of the n th electron and k(0,t), the counting 
process, is the number of electrons released during the time interval (0,t). Both 
k(0,t) and t n are random variables. A process like the one describe by (2.9) is 
called a shot noise process. In order to find the statistics of this process, 
specification of the statistics of the emission times {t n } and the counting 
process k(0,t) are required . 


15 


The electron count is the number of electrons that flow during a particular 
interval when detecting a radiation field with a photoresponsive surface. Given 
a constant irradiance field the probability of releasing k electrons from the 
photodetecting surface during the interval of time T is [Rl], 

P(k) = ~<T" k = 0 , 1 , 2 ,... (2.10) 


K has a Poisson distribution with parameter m. The parameter m is called the 
level of the probability and is directly proportional to the radiance level L that 
strikes the photodetector surface area A. 

m = 7L. (211) 


7 is a proportionality constant between the field intensity and the count 
intensity. It is important to mention that 7 is a function of the surface area A 
of the photodetector, the detector quantum efficiency, and the counting 
interval T. The characteristic function of the discrete random variable k is 
given by 

$(w) = E[eJ' wk ] = £ e> wk P k (k) (2.12) 

k=0 


00 

=e- m £ 


k=0 


m(e iw ) k 

k! 


exp(m(e^ w -l)J. 


It is well known that the characteristic function of a random variable can be 
used to find the moments of the corresponding random variable [Pi] , 
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E[k n j = -j n 


d n 


dw n 


(*k( w )] 


jw=0 


(2.13) 


From (2.12) and (2.13) it follows that the first two moments and the variance 
for the counting process k are 

E[k] = m (2.14) 

E[k 2 ] = m 2 + m 




= m. 


Rice [Rl] showed that the electron emission times in the photodetector are 
independent identically distributed random variables with uniform density 
functions over the interval (t,t+T), 


P(t n ) 


T t<t n <t+T 
0 elsewhere 


(2.15) 


Knowing the first order statistics of the counting process and the emission 
times, an expression for the density function of the Poisson shot noise process 
can be found, 

OO 

Px(x(t)) = TZ J <M w ,t)e“ jwx(t) dw, 

27r -00 


(2.16) 


17 


where the characteristic function of the shot noise, $ x (w,t), is given by 

*„<»,*) = Ele^M] = E k (E t Jexp(jw £ h(t-t„))| kj). (2.17) 

n=0 

Since the t„s are independent, 

= E k (nE, i |(ei” h < , -‘->]) = Ekl^fw.t)) 

n = I 

where <£ y (w,t) is defined as 

*,(".») = E,Jei” h <‘-*-»]. 

Taking the expected value with respect to k and making use of (2.12), 

00 1 

^xK*) = £ ^y(w.t)P(k) = exp([$ y (w,t)-l]m). 
k=o 

(2.17) is a general result given that a constant radiance is getting into the 
sensor. Without loss of generality, let us define the response function, h(t), as 
being equal to 

Sl 

T h 0<t<r h 

h(t) = 

0 elsewhere 


where q,, is electron charge and r h is the pulse duration time, then 
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4> (w)-l = ^r[exp(jw-^)-ll 
1 r h 


( 2 . 18 ) 


$x(w) = exp(nr h (e fk -1)) 


where, n = is known as the count intensity. For this particular realization 

of h(t), <l> x (w) is no longer a function of t. Again, the expected value of x and 
its variance can be found using the characteristic function <l> x (w), 


, _ . d$ x , _ _ m _ , , T 

E (x) - w=o - 9e n “ ~ 


( 2 . 19 ) 


2 _ n 9e _ , V7 


= 




=( vf ,L ' 


Observe, E(x) and <r x are both dependent on the received intensity. When the 
received intensity is high, the spread of the shot noise (its standard deviation) 
is small compare to E(x). In order to find an approximation to the shot noise 
density function let us define the normalized shot noise process as x n 


x n ~ 


_ x-E(x) 


( 2 . 20 ) 
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_.E(xi 

*x.( w ) = ^x(— )e a * 

U v 


(2.21) 


00 


1 ~a. 


= expf-jwy^ + E(i w ) q ( nr h) 2 )• 

q=l 


This last expression can easily be shown to converge to (2.22) as nr h —»oo 

,2 


\v 

KM - ex p(“yl- 


( 2 . 22 ) 


Since the normalize process x n is related to x by (2.20) it follows that x 
approximates a gaussian random variable with mean E(x) and variance 
provided nr h »l , 


Px( x ) * 




, -(x-E(x)) 

■ exp (. V A 11. 


2a‘ 




(2.23) 


Gagliardi [Gl], proved that this approximation is independent of the shape of 
h(t) as long as 

nr h »l (2.24) 


and 

OO 

J h q (t)dt<oo for all q > 2. 

-00 


nr h is physically related to the number of electrons emitted during the interval 
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(0,r h ) or in terms of the input radiance, 


n7 h 


= ( 



When the condition nr h »l is not satisfied, the discrete nature of the Poisson 
shot noise process will be evident at the photodetector output and a closed 
form representation for the density function (2.16) becomes very complicated. 

As was discussed in Chapter 1 the pixel radiance L from a given class and 
spectral band is a random variable. For a given pixel, L can be assumed 
constant because during the period of time that the measurement is done; there 
is no significant variation in the radiance being received by the photodetector. 
This suggest that in Landsat’s case what we have is a conditional Poisson shot 
noise process for every pixel observed on the earth surface. When this is the 
case, the characteristic function in (2.21) becomes a conditional characteristic 
function, conditioned over the random variable L. If the values assumed by L 
still satisfying the conditions of (2.24) then the conditional density function of 
the photodetector output will be (2.23) conditioned over L. For every observed 
pixel the output of the photodetector could be modeled as the pixel constant 
intensity level, the signal, plus a random variation from this intensity, the 
noise. This random variation or noise is normal with a zero mean value and a 
variance that is proportional to the (signal) intensity level. A model like this 
suggests looking at the shot noise process as a signal plus noise problem. There 
is one difficulty with this interpretation that must be kept in mind; it is the 
fact that both the signal and the noise are not strictly independent. 

It can be shown that the noise at the photodetector output is uncorrelated 
to the signal level; remember that two process can be uncorrelated but still not 
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independent. This is the case here. Since the noise in one channel is 
uncorrelated to the signal level in that channel, and different detectors are used 
for different wavelengths bands, the noise between the different photodetectors 
is uncorrelated. 

In this section the basis for the noise introduced by the photodetector was 
presented with an approximation for the density function of the photodetector 
output when the received signal level is high. 


2.4 Quantization Noise 


For Landsat satellites, the data undergoes an A/D conversion before it is 
transmitted to an earth ground station. A/D conversion involves the 
quantization of the analog signal to q=2 n levels , where n is the number of bits 
available. Basically, quantization is a mapping of the continuous data sample 
space into a finite number of selected values. Once the signal is quantized it 
can not be recovered exactly as it was before quantization. Since the 
quantization process introduces some fluctuation about the true signal value, 
this fluctuations can be regarded as an additive noise. The quantizer noise is 
defined as 

£ = x-x q , (2.25) 

x is the signal, and x q is the quantized signal. The mean square error 
introduced by the quantizer is given by [G2] 
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£ 2 = £/(x- Xi ) 2 p(x| i)dxP; = £/(x-Xi) 2 p(x)dx, 
i = ll, i = lli 


(2.26) 


where I ; is the i th interval range, p(x| i) is the probability density of the sample 
value when in I ; , and Xj is the quantization value of the i th interval. 

When a uniform quantizer is used and the conditional density function of 
the input is uniform over all intervals, the error e has a uniform distribution 
over one quantization interval A , 


P(*) = 


1 A A 

— -~ <e < — 
A 2 2 


elsewhere 


(2.27) 


and its mean square error is given by 


? = £i 
12 


(2.28) 


If the signal to be quantized has a normal distribution with zero mean and 
variance a 2 and a uniform quantizer with 2 q levels is used the mean square 
error is given by 


€ 


2 



(i + l)A 

/ (x-iA — — ) 2 e 2<73 d x + 

iA 2 


/ (x-^A + Y) 2 e 



.(2.29) 


Since the signal in Landsat satellites is usually modeled by a gaussian 
distribution, (2.29) is a more exact expression to describe the mean square error 
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introduced by the quantizer. From (2.29) it can be observed t 2 is a function of 
the signal variance. A computation of the mean square error (mse) and signal 
to noise ratio (SNR) was carried out for different standard deviations with the 
number of quantization levels fixed to 256, which corresponds to the Thematic 
Mapper quantizer system. This was done for gaussian and uniform density 
functions. The results obtained are shown in Fig. 2.3 and 2.4, which clearly 
show that for standard deviations in the range of .5 to 27 there is not 
significant difference between the mean square error introduced by quantization 
of a uniform or a gaussian density function. If we assume that the standard 
deviations of most of the input distributions to the quantizer in the TM system 
are in the range of .5 to 27, and the mean values are not near one of the edges 

A 2 

of the quantizer, then it is a good approximation to use as the mean 

12 

square error introduced by the quantization operation. 

By increasing the number of levels, making A smaller, we can make the 
mean square error as small as we want. In practice this can not be achieved 
because storage and bandwidth requirements impose an upper bound into the 
number of quantization levels actually used. 


log<mse) 
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Standard Deviation of Input Si anal 


Figure 2.3 Mean square error at the output of an eight bit 
uniform quantizer for input signals that are 
uniformly and gaussianly distributed. 


Signal to Noise Ratio <dB) 
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2.5 Johnson Noise 

Johnson noise , also known as thermal noise, is produced by the random 
motion of electrons in a conductive material. In 1924, Johnson performed a 
systematic observation of the thermal fluctuations of electricity in conductors. 
His observations were in complete agreement with the predictions done by the 
theory of thermodynamics and quantum mechanics. The theory predicts [G2] 
that the mean square voltage produced in a resistor R is given by 

(v-v m ) 2 = 4R(f)— ^ — Af, (2.30) 

e kT -l 


where: h = Planck’s constant, k = Boltzmann’s constant, v is the 
instantaneous voltage, T is the temperature in degrees Kelvin, Af is the 
bandwidth of the noise, and v m is the mean voltage. 

Thermal noise, black body radiation, and quantum noise can be developed 
from a few basic physical principles, as it was presented in a review paper by 
Oliver [Ol]. The voltage induced by thermal noise has been shown to have a 
normal distribution with zero mean and variance kTBR. 


1 { v* ^ 

^kTBR 6XPi 2kTBR ’ 


where B is the bandwidth of the noise . This result could be expected on the 
basis of the Central Limit Theorem, because the voltage generated in an 
electrical circuit by the thermal noise is the summation of individual short 
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current pulses produced by the electrons in the conductive material as they 
travel between collisions. For practical purposes it can be assumed that the 
thermal noise is white since the electron pulses produced by the collisions are 
extremely short in duration. Thermal noise is also observed in lossy dielectrics 
as the result of random thermal excitations of polarizable molecules, forming 
fluctuating dipoles. 
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CHAPTER 3 

MATHEMATICAL AND SOFTWARE MODEL 


3.1 Description of the Model 

In chapter two a statistical description of the noises in the Priority 1 list 
was presented. It is the purpose of this section to merge all of them together in 
a system that models the physical process of noise being added at different 
stages of the data acquisition process. In Fig. 3.1 a block diagram of the data 
acquisition system is shown, where all the noises are being referred to the 
radiance domain. This was done in order to avoid working with different types 
of units. We are not considering in this work the effects of the optics system in 
the signal. In vector equation form the model is described by 

Z = X(L) + U + V, (3.1) 

where, X, U, and V are independent random vectors. 

Z - is the sensor output, which is the signal to be transmitted to earth. 

L - is the radiance that gets into the sensor pupil after going through 
the atmosphere, 

L = r a L s + L p . 


X - is the photodetector output, it is a function of the signal level, 
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X(L,) = (r.L. + L p ) + N(0,KVr.L s +L p ). 


U - preamp noise, is modeled by a gaussian distribution with zero mean 
and diagonal covariance matrix. 

V - quantization noise, it is uniformly distributed with zero mean and 
A 2 

E v = 1, where A is the radiance range of one quantum level and I is 

12 

the identity matrix. 


The statistical properties of the input random process are modified at the 
different system stages. If we want to investigate the relationship between 
sensor parameters and overall performance we require a theoretical 
understanding of how the signal statistics are modified as it flows through the 
system. Although X(L) and V are not normally distributed, Z has been shown 
empirically [F3] to resemble a normal distribution. Arguments like the Central 
Limit Theorem are commonly used to substantiate this result. Since a 
multivariate gaussian distribution is completely specified by its covariance 
matrix and mean vector, we are interested to see how these random process 
parameters are affected by the system transformations. Another important 
descriptor to be considered, for each spectral band, is the correlation coefficient 
between the input and the output as a function of the system parameters. 

In Fig. 3.1 L is the actual radiance vector that enters the sensor pupil and 
X is the photodetector output vector in the radiance domain. As was shown in 
chapter 2 the density function of a single photodetector output can be 
approximated by (3.2) as the number of electrons emitted in the interval (0,r h ) 
becomes large, 
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P(Xi|Li) 


exp[ 


-(X-L,) 2 , 


V^l: 2ki 2 Li 


(3.2) 


kj is a parameter descriptive of the photodetector performance and it is related 
to the system parameters discussed in the shot noise section of chapter two by 
T 

k 2 = . The correlation coefficient between X ; and Lj is given by 

r h7 


Px^r 


E(X i L i )-E(L i )E(X i ) 


(3.3) 


where, 

E(X,) = ElJEx/X,! Lj)] = E(Lj) = l; ( 3 . 4 ) 

E(X; 2 ) = E L JE Xi (Xj 2 | L;)] = E(kj 2 L;+L 2 ) = kj^+L? (3.5) 

<4, = E(Xj 2 )-E(X;) 2 = kj 2 Li +*1 (3.6) 


00 

E(XjLj) = //XjLj P (Xj| Lj)p(Lj)dXjdLj = E(Lj 2 )=E?. 

-oo 


(3.7) 


So, Pxi,, is equal to 
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PX.L, “ 


*L. 


°L, 


°x. yVL, + cr l ' 


(3.8) 


This is a general result and does not rely on the asymptotic condition imposed 
over the shot noise process to obtain (3.2). Observe, p x l, is approximately 




equal to one if k; «-^rr ; this will correspond to an ideal photodetector. For 


a single spectral band the input-output correlation is 


^Z,L, ~ 


a Z?U 


(39) 


Since Xj,Uj, and V; are independent random variables 
<rz, = + *u. + *v, 


(3.10) 


E(z) = L; 


E(Z ; Lj) = E(Lj(Xi +U; +Vj)) = E(LjXi) = L; 2 


So, 




a z, 


(3. 11. a) 


Using the atmospheric model discussed in chapter two, Lj=7 a L s ;+L p , the 
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correlation coefficient between the input signal L si and Z; for a given spectral 
band is 


r,<7, 


^Z,L., - 




xAX. + i fci 2 (^ L 3i + L pl +4 +4, 


(3.11.b) 


As a test of the physical significance of (3.11. b) observe, for zero atmospheric 
transmittance Pzy. it = 0| and for no system noise, kj=<7u =<t v =0, the 
correlation coefficient is equal to one. Equation (3.11. b) clearly shows that the 
magnitude of the input-output correlation in a particular spectral band goes 
down as the average input radiance of the signal in that spectral band 
increases. 


For the definition of signal-to-noise ratio (SNR), using the variance of the 
signal, taken as t%o l m , over the variance of the noise, k^rJLsj+LpJ+a^+cryp 
we can express the SNR for each channel in terms of the input-output 
correlation as 


SNR 



(3.12) 


Recalling X i; U;, and V; are independent vectors we can state that the output 
covariance matrix is given by 

E z = £ x + E u + E v , (3.13) 


where the covariance matrix for U and V are diagonal matrices that are 
specified by system parameters. £ x is formed by diagonal elements <r x 2 , (3.6) , 
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and off diagonal elements P XiXj <x Xi <7 Xi where , 


PXJC, ~ 


ElXjXjl-Xj Xj 
a x, a x, 


Xi=L ; and ^ Xi 2 = k i 2 L i +<r L| 2 . 


(3.11) 


Using the model adopted in chapter two for a photodetector output, we have 


Xj(Lj) = L, + N(0,k iN /L[), 


(3.15) 


where the i th channel output is equal to the input radiance L ; plus a random 
variable with zero mean and variance kj 2 L ; . An expression for the cross spectral 
output correlation, Pxjt,’ can be obtained from, 

EtXiXj) =E[(L i + n i )(L j + n j )J =E(L i L j ) 


so, 


_ EfLjLjhL; Lj _ 


^L, 


(3.16) 


This is an interesting expression because it allows us to describe the 
photodetector output cross spectral correlation in terms of the input cross 
spectral correlation, ^ L(Lj . Combining (3.13) and (3.16) an expression for the 
output covariance matrix and mean vector can be obtained in terms of the 
different system noise sources, 
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E[Z] = E[L] = r a E[L s ] +L eq (I-r a ) 


(3.17) 


E„ = 


<r L i 2+k l 2L l +£r U 1 2+f7 'V 1 2 

PULfLfU ^L J 2 + k 2 L 2 + ^U 3 2 + ff V s 2 


(3.18) 


Ph x h t a h x °h % 


PLjLfLfU a h> + k n L„ + (r Vl 2 + a v . 2 1 


Since a covariance matrix is a symmetric matrix the upper diagonal elements 
in (3.18) were omitted. Observe that the off diagonal elements of E z are equal 
to the off diagonal elements of Ej. This implies a smaller correlation coefficient 
between spectral bands at the output of the system, because the expression 
inside the brackets in (3.19) is less than one. 

Pz,z, = Phi,, 


Vi, 


y / + k i 2 Ei +<7 ’u i +<T v 1 )( cr L 2 + kj 2 Lj +0u. +^Vj) 


(3.19) 


Summarizing, in this section an additive noise model has been adopted to 
describe the relationship between system parameters and the signal, and the 
output covariance matrix of this model was computed. It was shown that the 
higher the average radiance in a particular spectral band the lower the 
correlation between the input and the output of the system. Equations (3.11), 
(3.12), (3.17) and, (3.18) can be used in the evaluation or design phase of 
system parameters if statistical knowledge of the different kinds of signals to be 
observed is available. Having an analytical expression for the output 
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covariance matrices and the input-output correlation allows us to use 
parametric statistical pattern recognition techniques for the study of system 
parameters without having to actually generate the data sets representative of 
the noisy signal. 


3.2 Index of Performance 

If one is interested in the evaluation or design of system parameters a 
variable or variables against which any given parameter set can be compared is 
needed. Probability of correct classification or its complement the probability 
of error is a good choice for such an index of performance from a user point of 
view. Unfortunately, the Bayes minimum error can not be expressed in an 
analytical form because it involves a multidimensional integration over a 
complex boundary region, 

* = 1 ~ £ p K)/p( Z | W;)dZ. (3.20) 

»=i r, 


P(wj) is the a priori probability of class w; , p(Z| w ; ) is the density function of Z 
given that Z comes from class w i? and T; is the region in the multidimensional 
space where the Bayes criteria has classified the samples into class Wj. 

There are different algorithms which are commonly used for the error 
estimation. For example, random sampling techniques, such as Monte Carlo 
simulations, can be used for estimating the error. In this scheme, N random 
vectors are generated from the training samples statistics and they are 
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classified in accordance to the Bayes decision criteria. The number of 
incorrectly classified samples divided by i\ gives (he maximum likelihood 
estimate of the Bayes error. This estimator is known to be unbiased [Fl], The 
error estimation to be obtained from a Monte Carlo simulation is a random 
variable which is dependent on the number of samples used. For accurate 
estimates a large number of samples are required which implies large amounts 
of CPU time. An algorithm was presented by Fukunaga and Krile [F2], for a 
two class problem, which transforms the multidimensional integration of (3.20) 
into a one dimensional integration in the frequency domain. This algorithm, in 
theory, is exact. Although, small errors may be introduced by its numerical 
implementation, compared to Monte Carlo simulations, it is more accurate, and 
more efficient in machine time. 

Statistical separability can be used as an alternative index of performance 
when the probability of error results are inconvenient to compute. A well 
known statistical distance function between two gaussianly distributed classes 
is the Bhattacharyya distance[BlJ, 


V = 


1 _ E,+E 2 1 I 2* E ‘ +E ^I 


y 2 y 


(3.21) 


which has been commonly used for feature extraction applications [S3]. It is 
appealing to use the Bhattacharyya distance as an aid in the design or 
evaluation of system parameter. Although, there is not an exact relationship 
between probability of error and the Bhattacharyya distance, an upper and a 
lower bound exists for the Bayes error in terms of fi [Fl], 
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|[l-v/l-4£ u 2 ] <£<<„ = je-“. 


(3.22) 


Based on simulation results Whitsitt [Wlj suggested the usage of the 
complementary error function of s/fy 7 as an approximation of the Bayes error. 
When the two classes under consideration have equal covariance matrices and 
equal a priori probabilities, the approximation turns into an equality. If this is 
a good approximation it has to be bounded between the upper and the lower 
bound of the Bayes error, given by (3.22), for all values of /i. A proof of the 
consistency of this approximation follows. 


(3.23) 


The complementary error function is defined as 


OO 


-x s 


Q(x) = (tit 2 dx - 


It is known Q(x) has an upper bound given by [VI], 
_xi 

Q(x) < — e 2 x > 0. 

2 


(3.24) 


Replacing x by >/2fi in (3.24) it is easily proved Q(v/2p) is less than the upper 
bound given by (3.22). The lower bound consistency can be proven by putting 
the left hand side of (3.22) in the form 

1 < 2Q(\/2 Ji) + \f I~e -2/4 = g (/*). (3.25) 


Observe that at p equal zero or infinity the equality is satisfied. Since g(//) has 
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these boundary conditions, if it ever goes below the value one, it will have a 
minimum point. Taking the derivative of g(ji), using the Leibnitz rule of 
derivation under the integral sign, and setting it equal to zero, we obtain 

-l -l 

g'(/z) = 2 +e-^(l-e- 2 ^) 2 = 0. 

V27T 

Since we are interested in the relative minimum in the open interval (0,oo) the 
previous expression can be simplified into 
e- 2 "[7r// + l]=l, 


which can be solved interactively obtaining the value /i — .4221 . Substituing 
this value in (3.25) we obtain 
1 < g(.4221) = 1.113 . 

This proves the approximation is consistent with the bounds established on 
the Bayes error by (3.22). The advantage of using Q as an approximation 
of the Bayes error is that it does not require extensive computer computations. 
In Fig. 3.2 the relationship between Q(v/2 jz) and the upper and lower bounds 
on the error probability is shown. 


Probability of Error 
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Figure 3.2 C 
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3.3 Software System 

The next step was to construct a software system which could adequately 
simulate noise from each of the types of sources to be considered. In order to 
simplify the problem, only two classes with equal a priori probabilities were 
simultaneously considered, and the output distribution of these was assume to 
be gaussian. Since real data is likely to have characteristics that are highly 
dependent on outside and uncontrollable elements, such as deviations from the 
gaussian assumption, it was decided to use the covariance matrices obtained 
from training samples as the input signal. Thus a system of computer 
programs was written to allow for the use of a LARSYS data statistics file as 
input and to use it to simulate an unadulterated signal source. In general the 
required inputs to the software system are: 

1- Atmospheric optical thickness vector and solar zenith angle 

2- Duntley’s equivalent radiance vector, L eq 

3- k vector, with shot noise system parameters 

4- Covariance matrices for preamp and quantization noise 

5- Spectral covariance matrices and mean vectors for the two classes to 
be considered 

6- Weighting factors vector for the selective adjustment of the different 
noise sources 

The software output consists of: 

1- Bayes minimum error for each output class and overall 

2- Bhattacharyya distance between two output classes 

3- Complementary error function of \/2/7 

For computing the pairwise error, the algorithm developed by Fukunaga and 
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Krile was used. A list of the software programs can be found in Appendix C. 


3.4 Selection of the Data Set 

As an initial data set airborne multispectral scanner data of Segment 210 
of the 1971 Corn Blight Watch Experiment was selected. The scanner 
instrument used to collect this data was the Michigan M-7 airborne scanner. 
The selection of this data set was influenced by the following factors: 

1- Ground truth for every pixel should be available. This is important 
both for deriving good quality training samples and for accurate 
determination of performance. 

2- The characteristics for this data set, e.g. the number of spectral bands 
and the ground spatial resolution, are representative of possible future 
spaceborne scanners. 

3- The informational classes in the ground scene are adequately complex 
to be representative of possible future applications and to provide a 
significant challenge against which to test future scanner system designs. 
Out of this data set we used a subset of classes that gave good 
classification performance when no noise was present, but was sensitive 
to the noise levels introduced by the sensor system. The criteria 
adopted to choose these classes was based on the coincidental spectral 
plot from the training samples of segment 210 . 

It was decided to use bands and noise levels which approximate the 
Thematic Mapper in order to enhance the practical benefits of the study. 
Channels 1,5,7,8,10 from segment 210 data were the ones which had the best 
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match with channels 1,2,3, 4, 5 in the Thematic Mapper. The Thematic Mapper 
noise equivalent radiance curves for quantization, preamp, and shot noise were 
obtained from NASA’s Goddard Space Center for the first five spectral bands 
in the TM, see Figure 3.3 to 3.7. This curves clearly show the dependence 
between shot noise and the signal level. 


3.5 Simulations Performed and Results 

In this section it is shown how, by using the model developed, we can 
investigate the relationship between classification accuracy, sensor parameters, 
class statistics and the atmospheric effects. In particular, we want to study if 
by slightly adjusting the noise levels introduced by the Thematic Mapper 
sensor system we can make a significant change in the recognition accuracy. 
Another parameter to be modified is the atmospheric visibility. 

A series of different computer simulations were conducted using the 
developed software system. In these runs the data set used came from the Corn 
Blight Watch experiment. We considered this data set as representative of a 
noise free data set and for the purpose of the simulation we assumed it was 
collected under perfect atmospheric conditions by an ideal TM kind of 
instrument. Using the software system we can corrupt the signal statistics, 
based on the chosen system parameters, and quantify the degradation in terms 
of the index of performance proposed in the previous sections. 

In the first experimental run the inputs were: 

1- Equivalent _ radiance vector, L eq = 150(1, 1,1, 1,1) T , the particular 
choice of this vector value was arbitrary. 


Noise Equivalent Radiance 10 3 W/cm d -sr 
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nl — shot noise 
n2 = pream noise 



p 

Signal Radiance mW/cm^-sr 


Figure 3.3 Thematic Mapper noise equivalent radiance in the 
•45 - .52 um band. Channel gain is set to have 
1.06 raW/cm2-sr = 256 bits. 


Noise Equivalent Radiance 10""^W/cm^-sr 
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nl = shot noise 
n2 = pream noise 
n3 = quantization i 



p 

Signal Radiance mW/cm -sr 


Figure 3.4 Thematic M aDDe r noise equivalent radiance in the 
.52 - .60 urn band. Channel gain is set to have 
2.54 mW/cm^-sr = 256 bits. 





Noise Equivalent Radiance 10 5 W/cm 2 
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nl = shot noise 
n2 = pream noise 
n3 = quantization noise 



p 

Signal Radiance mW/cm -sr 


Figure 3.6 Thematic Mapper noise equivalent radiance in the 
.75 - .91 urn band. Channel gain is set to have 
3.26 mW/cm^-sr = 256 bits. 


Noise Equivalent Radiance 10 °W/cm 
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d 1 — shot noise 
n2 = pream noise 



figure 3.7 Thematic Mapper noise equivalent radiance in the 
1.55 - 1.75 um band. Channel gain is set to have 
.64 mW/cm^-sr = 256 bits. 
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2- Atmospheric optical thickness matrix for a haze atmosphere. Eight 
different metereological ranges were considered, see chapter two. 

3- Thematic Mapper noise equivalent radiance system parameters for 
shot, preamp and quantization noise; this was obtained from NASA. 

4- Covariance matrices and mean vectors for the two input classes, see 
Fig. 3.8. This matrices correspond to the subclasses soyl and soy2 
obtained from the training samples selected by LARS staff. 

In Fig. 3.9 is shown the results of how the Bayes error changes for 
different atmospheric visibilities, which correspond to a haze regime. Observe 
that for high visibilities the error due to preamp noise alone is significantly 
larger than that due to shot noise and quantization noise. This is in complete 
agreement with the kind of theoretical result expected from (3.17) because the 
atmosphere attenuates the power in the signal and introduces a path radiance 
component into the signal. This causes the two input signals to approach one 
another making them more susceptible to system noise. Fig 3.10 shows that the 
same relative results would be obtained when Q(\/2/i) is considered, although 
the values of Q (\/2ju) are always least when compared to the values of the 
error. Figure 3.10 is a smoother version of Fig. 3.9; this is explained by the 
fact that the function Q(V2ji) is an approximation of the actual error and can 
not track perfectly the fluctuations in it. Quantization noise did not show a 
significant variation in the recognition accuracy for the range of atmospheres 
considered and was the less degradating noise source in the system. 

The results in Fig. 3.11 correspond to an atmosphere with a horizontal 
visibility of eight kilometers. Here it is shown how the Bayes error changes 
when the shot noise parameters, K, are adjusted by a variable c from zero 
times its original or actual value to twice its actual value. While this 
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M t = [120.98, 114.74, 81.96, 144.39, 162.24] T 


Ei 


143.06 

142.52 167.5 
179.29 209.13 276.90 
-19.62 -34.3 -51.5 66.5 
81.7 88.01 112.95 1.96 73.81 


M 2 = [115.36, 109.26, 76.94, 141.65, 149.94] t 


e 2 


16,8 

5.11 13.94 

2.61 6.11 7.72 

-.99 -3.06 -4.83 25.5 
-1.35 -7.11 -2.90 7.43 39.92 


Figure 3.8 Statistics for subclasses soy 1 and soy 2. 


parameter was changed the rest of the system parameters were set to their 
original values. There was a change of approximately five percent between the 
maximum and the minimum of the Bayes error in Fig. 3.11. Figure 3.12 shows 
the same type of result when the preamp noise level is the adjusted parameter. 

In Figure 3.13, the adjusted system parameter was the number of bits 
used by the quantizer for representing the dynamical range of the incoming 
signals. When this result is compared to Fig. 3.9 it is observed that if five bits 


Probability of Error </;> 
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Met ereo logical Range 


Figure 3.9 Probability of error versus metereological range. 
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Metereological Range Ckm> 


Figure 3.10 Bayes error approximating function versus mete- 
reological range. 
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were used, the quantization noise turns out to be the most predominant source 
of noise. For the particular set of classes considered in Fig. 3.13, an 
improvement in the number of bits from eight shows no significant difference in 
the probability of error. 

The purpose of this section was to show how the software system could be 
used to evaluate the response of a multispectral scanner under different 
operational environments and different sensor parameters. The results obtained 
here were specific for the two input classes selected. 


3.6 Removal of Atmospheric Effects and Noise 

As was observed with regard to Fig. 3.9, th? atmosphere has an important 
role in the degradation of class recognition. Based on the noise free atmospheric 
model discussed in Chapter Two, it is reasonable to ask; if the atmospheric 
transmittance and path radiance are available to us will it be possible to do an 
atmospheric correction on the data such that the classification results can be 
improved? In order to answer this question, lets first look at what an 
atmospheric correction would consist of. To correct for the atmosphere an 
inverse type of operation is needed. This will involve changing the data from 
the digitized domain into the radiance domain, subtracting from it the 
estimated path radiance and then multiplying by the estimated inverse 
atmospheric transmittance matrix. By looking at (3.1) we can easily see that 
the resulting expression will be 
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Y = r a -‘(Z-L p ) = L,+t^N(L). 


(3.26) 


where N represents all the noises introduced by the sensor system. For data 
taken under good atmospheric conditions the second term could be negligible 
compared to the signal term, L s . Now, if the classifier used is a Bayes 
minimum error pixel classifier the answer to the question is no. This is because 
as was shown in Chapter Two, this type of operations does not affect the 
percent of overlapping between the possible classes. 

If it is desirable to estimate the reflectance of the ground based on the 
available data, the digitized data will have to be converted to the radiance 
domain. Once in the radiance domain all atmospheric effects present in the 
data have to be removed as mention earlier and then multiplied by the inverse 
solar irradiance matrix measured at the earth surface. Much has been done 
with this kind of calibration procedure [S2], which sounds very reasonable if a 
photo-interpretation method is to be used. But again, if the main concern of 
the user is to do an accurate estimation of the classification based only on the 
frame of data available, and using a quantitative approach then it makes no 
difference in the classification results to correct for the atmosphere or go to the 
reflectance domain. 

Maxwell [Ml] has made use of spatial filters as an alternative method of 
reducing the spectral noise and in that way improving the classification 
accuracy. In particular he made use of a moving average window which 
showed a classification improvement from 18.8 percent of error to 3.6 when it 
was applied to a noisy data set. He made the observation that although there 
was a significant improvement in the classification accuracy the edges of the 
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final images were blurred by the filter. Since a linear spatial filter like the 
moving average window will tend to blur the edges in the image, defined at the 
boundaries of one class to another (introducing possible misclassification), we 
decided to seek another kind of filter which could preserve a good classification 
accuracy in the overall image frame. A non-linear filter known as the median 
filter seemed to have these attributes [H2]. 

The implementation of a median filter can be understood in a similar way 
as the moving average window, but in the median filter case the pixel falling in 
the center of the window is replaced by the median value obtained from all the 
pixel values falling inside the window. Some preliminary results were obtained 
using a recursive separable median filter with a window of width three. Table 
3.1 summarizes the results obtained when a four band multispectral image was 
corrupted with additive gaussian noise with standard deviation equal to seven. 
This noisy image was filtered with a recursive separable median filter and then 
a classification was performed and compared to the classifications obtained for 
the original image and the noisy image. The classification improvement after 
filtering was greater than twenty percent overall. Figures 3.14 - 3.16 show the 
classification results. One of the disadvantages of this filter is that the output 
statistics are far from trivial to obtain and closed form expressions are very 
difficult to find. Future work is planned to be done with this kind of filter and 
to compare it to the moving average window and other proposed non-linear 
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Table 3.1 Classification results for a separable recursive 
median filter. 


Class 

Percent of Correct Classification 

Original 

Noisy 

Filtered 

Soybeans 

96.9 

42.0 

73.2 

Corn 

97.0 

62.3 

86.2 

Oats 

95.3 

51.9 

74.9 

Wheat 

97.3 

63.3 

83.4 

Red cl 

92.6 

51.3 

76.0 

Alfalfa 

94.9 

79.3 

88.0 

Rye 1 

95.7 

50.5 

86.2 


* 



Figure 3.14 Classification results of the original data. 
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Figure 3.15 Classification results of the original data 
plus white additive gaussian noise, N(0,7). 
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filters. 

3.7 Observations 

In this work a model for the data acquisition system in a multispectral 
scanner system like the one utilized by the LANDSAT satellites was presented. 
Since the shot noise introduced by the photodetectors in the sensor system is 
signal level dependent, an atmospheric model was adopted which could 
adequately describe the amount of radiation that gets into the sensors based 
on the atmospheric transmittance. An analysis was carried out to find the 
output spectral statistics in terms of the input signal statistics and the system 
parameters. This was integrated into a set of fortran programs that can be 
used to estimate the classification performance when supplied with the class 
statistics, noise levels introduced by the sensor system, the atmospheric 
transmittance, and the atmospheric path radiance. 

F urther topics to be considered in the improvement of this model are: 

1- Relaxation of the assumption of relatively high input radiance 
arriving at the photodetector input, which allowed approximating the 
photodetector output by a gaussian distribution with parameters being 
signal level dependent. 

2- Implementation of a non-constant atmospheric model. 

3- Use of a path radiance model that considers the cover class under 
observation. 

A library of programs called the Unified Scanner Analysis Package (USAP) 
had been developed at the Laboratory of Application of Remote Sensing 
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(LARS) to simulate a multispectral scanner in a more global way than the one 
presented in this thesis. This package has a weak point because it does not 
consider any atmospheric effects and it assumes the sensor noise to be the same 
in all spectral bands. The work done in this study could be easily appended to 
the USAP system developed at LARS in order to obtain a better simulation of 
the multispectral scanner system. 
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APPENDIX A 
THEMATIC MAPPER 


The Thematic Mapper (TM) is a scanning optical sensor aboard the 
Landsat 4 satellite, which is in a near polar sun- synchronous orbit at an 
altitude of 705 kilometers, with an earth swath scan of 185 Km wide. The TM 
is a second generation of earth resources space sensor. It works on the same 
basic principles of the Multispectral Scanners (MSS) aboard Landsat 1, 2, and 
3. It is an image forming system where each pixel in the image is a vector 
consisting of a set of measurements from selected wavelengths regions in the 
spectrum. The TM was design so that it will achieve higher imagery 
resolution, sharper spectral separation, improved geometric fidelity and greater 
radiometric accuracy and resolution, compared to previous scanners. The 
spectral bands used by the TM are 


band 

range fim 

description 

1 

.45-. 52 

blue 

2 

.52-. 60 

green 

3 

.63-. 69 

red 

4 

.76-. 90 

near infrared 

5 

1.55-1.75 

middle infrared 

6 

2.08-2.35 

middle infrared 

7 

10.4-122.5 

thermal infrared 
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These bands were selected because they have proven to be helpful in the 
prediction of vegetation mapping, which is the principal mission of the 
Thematic Mapper. Silicon detectors are being used for the first 4 bands. The 
two middle infrared and thermal infrared bands use indium antimonide and 
mercury cadmiun-telluride detectors, respectively. The detectors of the latter 
three bands must be at low temperatures in order to preserve high SNR. 

Harnage and Landgrebe [Hi] have stated that classification accuracies are 
acceptable for most remote sensing applications when the fields of the scene are 
greater than 8 Instantaneous Fields of Views (EFOV) in size. With the TM 
IFOV of 30 meters, adequate crop estimation is more approachable in countries 
where the field sizes are smaller that the ones in USA. 

A nominal equatorial crossing time of 9:30 AM has been chosen for 
Landsat 4 because it will provide less cloud cover and topographic information 
can be obtained from the shadows. It follows that a nominal equatorial solar 
zenith angle of 37.5 degrees will be achieved by the Thematic Mapper. For a 
detailed description of the Thematic Mapper design predicted performances see 
[HI, SI], 
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APPENDIX B 

NON CONSTANT ATMOSPHERE 

The non constant atmosphere considered here is one where the 
atmospheric transmittance and path radiance varies along the satellite path. 
This could be induced by inhomogeneity of aerosol and water vapor 
concentrations in the atmosphere. A model presented by Kiang [KlJ that 
considers this fluctuations is given by 

L = r a (t)L s + L p (t) (1) 

where r a and L p are the transmittance matrix and path radiance vector, 
respectively. Both of them are considered to be random and are functions of 
the optical thickness fc which is also a random vector with mean value t 0 . Let 
r 0 = r a (t 0 ) and L po = L p (t 0 ). Expanding L in a Taylor series of t about t 0 . 

L = r 0 L, + L p- + |(^) 0 L, + (^) 0 ](t-to) + H.O.T. (3) 

Assuming that the fluctuations of the optical thickness are small we can neglect 


the higher order terms and make the following approximation, 
L - r o L s + L Po + 7 


( 4 ) 
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which resembles the constant atmospheric model presented in Chapter 2 with 
the exception that there is an additive term 7 which is due to the fluctuations 
of optical thickness. This term can be considered as a first order 
approximation to be the atmospheric noise. It has an expected value of 

EM = 0 (5) 

It then follows 

E|L] * r„E|LJ +L Po 

El =« r 0 E u r' + E, 

This model was not used in the present work because the statistics of the 
random variable 7 were unknown. 
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APPENDIX C 

SOFTWARE LISTING 


program mss 

* This program simulates Ithe a Latin tie a] transformations that 

* two Input signals undergo at different stag as of a 

* multlepectral acanner system The signals for each data 

* are characterised by the covariance metrical and mean vectors 

* The program eatl mates the probability of mis class location and 

* the Bhattacharyya distance between classes The Inputs to the 

* program are 

* lname=output file name 

* tao=otmospherlc optical thickness values 

* a ^values that will adjust the actual nolle levels. 

* mr=metereologlcal range 

* factor- increment used by numerical Integration 

* dimensionality to be used 

* ko=ectual dimensionality of input statistics 
9 Index =lth coefficient In e to be modified 

* value l=atarUng value of a(lndex) 

* valueZ=ftnal value of a(lndex) 

* step 1=# of intervals used 

9 cnanae(iat2)= input file with signal statistics 

* Libraries used (IMSL) are used 

WW tWIM I MM M W M t W 

/ > character*6 cname(Z).fname 
‘ real ks(25.2).mz(25.E).kys(21.Z).eys(6.2).lp(8).t(S) 
real s(25.4).rmax(6).rmln(6).a(3).tao(5).kl(6).k2(6) 

data kl/2 308e-4.2 1533e-4.2.363e-4.1 76e9e-4.1.757Be-4.n.7B7Be-4/ 

data k2/ 74e-b..53e-5. 4be-b..3ft«-b. 72e-b..»2e-5/ 

data lp/150 .150 .150 . 150 .150 .150./ 

data rmax/ 1.05e-3.B &4e-3.1 46e-3.3.26e-3..54e-3..4B7e-3/ 

data rnftln/0 .0..0 .0 .0 .0 / 

data cname/^claaBlVclassZV 

* read control flle 
read(5. 9 ) fname 

crpen(E.m« -fnamn) ■ 

rewind 2 

1000 read(5.*.end=300) tao 
read (5, •) a.mr.factor.k.ko 
re ad(5.*) Index, value l.volue2,stepl 
step=(value2-value 1 ) /slap 1 
print* ("fee tor=".fB 4)'. factor 
wrlte(6,*) 
kd=ko*(ko+l)/2 
kdl=k»(k-M)/2 

* read Input statistics for two classes 
do 1 1=1,2 

open( l.JUe=cname(l)) 
rewind 1 

read(l, # ) <kysUl).J=I.kd).(eysUU.)=l.)«>) 
close(l) 

1 continue 

* set values of weighting eoeff 
elude x-valuel 
lfltep=nlnt(stepl) 

if (value. 1 eq value 2) istep=0 

do 2 iqv=0.1fltep 

a(lndex)=alndex+ lqv*step 

wrlte(6,6) Index, value l,value2,stepl.a(lndex) 

8 format(13,4f8 2) 
wrlte(fl,*) 


* form cov matrices for shot, pres m quant noise 

the ta=(37. 5/350 )»2 *3 14158 

do 3 1=1 Jt 

la=l*(l-l)/2 4 1 

■( la. 1 ) =a ( 1 ) *k 1 (1) *s qrt(255 . / (rm ax(l>-rmln(i) )) 
s(la,2)=a(2)*sqrt( (e55 /rmax(l)) *k2(l) ) 

*(la.3)=(a(3)**2)/12 
t(l)=erp(-l *tao(l) /cos (theta)) 

3 continue 
lf(lqv eq.0) then 

print* ("mete re ologlcal range : **.12)*. mr 

wrlte(6.*) 

print* ("atmospheric transmittance **.10f6.3)\(t(l)J=lJc) 

write (6.*) 

en dlf 

do 500 le=l,2 

* farm covariance matrix for ■'.fin a. 

do 1001x=l.k 
pathr= lp(lx) *(l-t(lx)) 
mi(U.lc)=t(u)*eTi(U.lc) * pothr 

do T5ly=ljx 
klr=lx*(U-l)/2 4 ly 
U(U oq.ly) than 
r=t(U)«ky»(kk.ic)n(ly) 

k*(kk.le)=« 

else 

ki(kk.le)=Ua)»ky»(kk.lc)*t(ly) 

endlf 

75 continue 
100 continue 
500 continue 

c call utwsm(* , ilgJ’\4.kt(l.l)Jt.l) 
c call u»wsm( ,, *lg2'’.4.kf (1,2) Jc.l) 

call dlag(kt(l.l).kt( 1 . 2 ).mE(l.l).mE( 1 . 2 ).kitdUactor.a(te*lfi)) 
2 continue 

go to 1000 

300 print'C'end of flle'O' 
stop 
end 
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subroutine dlag(sigl,sig2.mvl,mvg.nd.ndl.fBCtcr.vlndcw) 


• This subroutine does e simultaneous dlag email cation on the two 

• Input c overlance matrices The transformation used Is also 

• applied to the mean vectors It computes the Bayes error. 

• computes the Bhattecharyye 

• distance and the upper bound of the Bayes error etabllshed by 

• It. 

• SUBROUTINES CALLED. 

• subroutine error end subroutines defined la the MSL library 

reel mvi(nd),mv2(nd).ilgl(iidl).slg2<ndl) 
reel 1(10. 1 0), tt( 10.10). irk 1 (200) 
real phl(10.10).wk(2n0).wkmO0.10).eval(50),dlit(50) 
real evec(10.10).wke(10.10),res(10) 

)obn=2 

lr=10 

ao 12 lx=l.nd 
do 12 ly=l,nd 
12 wkm(lx.ly)=0 0 

• flnd eigenvectors and eigenvalues of elgmal -> elgenv.phl 
call aigrs(sigl.nd,)obn.oval,evee.U,wk.l«r) 

lf(wk(l) gt 100 ) then 

prlnt’C performance lndexl=”.fl4.5)*,wk(l) 
stop 
endlf 

• form eval**(-.5) 
do 30 li- l.nd 

30 wkm(lx.lx)=eval(lx) ## (-.5) 
c call u»wfm("evscl".5.evwc.l0.nd.Bd,2) 
c call uswlm("evnl**-.5".ll.wkm.l0.iid.mL2) 

• farm T= eva)*«(-.5) • erec-t 
do 40 lx=l.nd 

do 35 iy=ljid 
sum=0 0 
do 33 1=1 .nd 

■um=sum+wkm(lx,l) •evee(iy.l) 

33 continue 
t(lx.ly)=sum 
U(ly.lx)=sum 
35 continue 
40 continue 

c call uewfm(T\i.ui0.mLiuL2) 

• form t*slg21t -> tt 

call vmulsf(Big2.nd.tt,nd.lO,wkm.lO) 
call vmulffO-.wkm.nd.zuLnd.lO.lO.tt.lO.ler) 
c call uswtm(**k=T*slg2*TT*\ll.tH0.nd.iid.E) 

call vcvtfs{tUnd. lO.wk) 
c call uewsmC'K’.l.wk.nd.B) 

• And eigenvectors and elgevalues of wk 
call elgrs(wk.nd, John, aval evec.ic.wkl.ler) 
lf(wk(l).gt 100) then 

prtnt*( "performance lnde*2=‘.fl4 5)\wkl(l) 
stop 
endlf 

c cell uswfm("evee of k*’.B.evec.l0.ndmd.2) 
c prlnt’C eigenvalues ofK".10fl4 5)\(eval(l).l=l jsd) 

• form phl=evee-t • t 

call vmulfm(evec.UmLnd.nd. 10. 10. phi. 10. lor) 
c call uswfm( "matrix trans ".13. phi lO.nd.nd.E) 

• form dlst vector 
dc 60 1=1. nd 

mv2(l)= mv2(l)-mvl (1) 


rum =0.0 
do 50 )=l.nd 

tt(U)=phiUl) 

50 continue 
60 continue 

call vmulff(t,mv2.iuLnd.l. 10.50. dist. 50. ier) 
c call uswsm ("slg 3 ".4, clg l.nd. 1) 

c cell uswsm("slg2~.4.slgE.nd.l) 

• form Identity 

c cell vmulaf(slgl jid.tUnd.10.wkin.10) 
c eall vmulff(uwkm.nd juLnd. 10.10. wka.10.ler) 
c call uswlm( "ldenUty*’\9.wke.l0.nd.nd.l) 

* form eigenvalue matrix 

c call vmulsf(slg2jid.tUnd.l0.wkm.l0) 
c call vmulfT( t.wkmmtLnd.nd. 10.10. wka.10.ler) 
c call uswfm("elgenvaines*".l2,wka. lO.wLnd.l) 


* compute Bhattacharrye distance 
sdlst=0 

do 100 l=l,nd 

dd=dlst(l)**2/(E*(l ♦eval(l))) 
dd=dd* log((l.+eval(l))/2 ) -.5*log(evatfl)) 
sdlsteedlst+dd 
100 continue 
edlst=.5*Bdl>l 
el=txp(-l ••diet) 
res(l)=sdlst 
res(2)=e 1*100 
res(3)=el*50 

res(4)=erfe( )*50. 

c prtnt’C'u = ••.16.3)'.«dl» t 
e el=exp(‘l.*rdiet) 

c print*' l < .U V10x."e <*\fB 3)'.al*100..100 *el/E. 

• com?L. 'hciuaT error 

call emr(s^USlft.tUdlsUevaLfactorjid.al > e2.a.nf) 

re*(B)--*1 

re»(fl)=e2 

res(7)=e 

res(B)=vlndex 

res(B)=Hoat(nf) 

wrlle(6, 1000) (res(l)J=1.0) 

1000 format( "results ~.5fl2.3) 
write (B.*) 

wrtte(2,*) (rws(l)J=l.Q) 

return 

end 
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tub routine error<a.b.h.d,elg.faetor.iHilm,el.a2,e,l£) 

• Thi* eubroutlne 1» called by “ding"- It mUbbUi the 

• probability of error between two cla*ee* by performing 

• a numerical Integration over a single variable In the 

• frequency domain, tea Fukunaga and Krlie. 

real a(ndlm.2).b(ndim,2).h(2idlm.2) l d(ndim).elg(xidlin) 
c print* (“enter factor ndim elg d ")* 
c read(5.*) factor.ndlm.elg.d 
wrlte(B,*) 

print* (“eigenvalue* ".10fB.3)\(elg(l).l=l.ndim) 
write (B,*) 

print* (“dlitance vector*", 1 Of B.3)’,(d(l),l=l.ndlm) 
wrlte(6.*) 

11 lormat(10fB.3) 
lc=0 
el=0 
e2=0 

pl=3 1415B2B54 
vl = l 
v2=.25 
v3=.5 

do 5 1=1, ndim 
dl*t=d(l) 

a(U)=l -l-/elg(l) 
a(L2)=elg(l)-l 
b(U)=(l /#lg(l))*dl»t 
b(1.2)=»qrt(elg(l))*dl*t 
bl2=b(l.l)**2 
b22=b(l,2)**2 
lambda=log(elg(i)) 
h(tl)=bl2/(l -a(l,l)) + lambda 
h(1.2)=-l *b22/(l +a(LZ)) + lambda 
5 continue 

1 continue 
lc-=lc+l 
w=lc*f actor 
w2=w"2 
pmag 1=1. 
pang 1=0. 
p mag 2=1 
pang 2=0 

do 10 1=1. ndim 
aie=a(l,l)**2 
a22=a(l,E)**2 
blB=b(U)**2 
b22=b(i.2)**2 
awl=a(l.l)*w 
aw2=a(l,B)*w 
alw2=a(l.l)*w2 
a2w2=a(L2)*w2 
dml=vl4-w2*alE 
dm2=vl+w2*n22 

duml= atan(awl>-w*((h(i.i)+bl2*Blw2/dml)) 
dum2= atan(aw2)-w^((h(1.2)+b22 •oEwE/dmE)) 
pang 1 =pang 1 ♦ duml 
pang 2= pang 2 -f dumE 

pmag 1 =pmag 1 *(dml ••(-v2) •e*p(-v3*bl2*w2/dml)) 
pmeg2=pmag2*(dm2**(-v2)*azp(-v3*bE2*w2/dznE)) 
10 continue 
pangl=.B*pang 1 
pang2=.9*pang2 
phll =(pmag 1 / w) *eln(pang 1) 
phl2=(pmag2/w)*uln(poagg) 
c print '("ang, mag.phie,w”,4f 14. 5)‘.pang2,pmag2,phl2,w 
el=el+phll 
e2=e2+phl2 
pmag=pmagl+pmag2 
If (pmag gt..0001) go to 1 

el= .5 - (l./pl)*el •factor 

el = (1 - el)*100 

e2= ( 5 - (1 /pl)*e2«faetor)M00 

e= (el+e2)/E 

e prlnt*("el.e2.e.n".3l7.2.lfl) , ,el,e2,e,lc 
return 
end 


